{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time domain reflectometry, measurement vs simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrate the use of frequency to time domain transformation by comparing measurements and simulations of a [microstripline](#Microstripline) and a [microstripline with stepped impedance sections](#Stepped-impedance-microstripline).\n",
    "\n",
    "The simulation data is generated by `skrf` using a simple transmission line model for connectors and each impedance section. To achieve a reasonable agreement between measured and simulated data, the dielectric permittivity as well as the connector impedance and delay are extracted by optimization. The code for the [simulation](#Simulation) and [parameters optimization](#Parameters-optimization) is given at the end of the example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from IPython.display import Image\n",
    "from numpy import absolute, log, sum\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "import skrf\n",
    "from skrf.media import DefinedAEpTandZ0, MLine\n",
    "\n",
    "skrf.stylely()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data into skrf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The measurement was performed the 19th February 2018 on an Anritsu MS46524B 20GHz Vector Network Analyser.\n",
    "The setup is a linear frequency sweep from 1MHz to 10GHz with 1MHz step, 10000 points.\n",
    "If bandwidth 1kHz, 0dBm output power, no smoothing, no averaging. Full two-port calibration with eCal kit.\n",
    "\n",
    "Considerations about [time resolution](#Time-resolution) and [range limitation to avoid alias response](#Measurement-range-limitation-to-avoid-alias-response) are at the end of this document."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load all measurement and simulation data into a dictionary\n",
    "meas = skrf.read_all_networks('tdr_measurement_vs_simulation/measurement/')\n",
    "simu = skrf.read_all_networks('tdr_measurement_vs_simulation/simulation/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DC point extrapolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The measured and simulated data are available from 1 Mhz to 10 GHz with 1 MHz step (harmonic sweep), thus we need to extrapolate the dc point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['P1-MSL_Stepped_140-P2',\n",
    "         'P1-MSL_Thru_100-P2']\n",
    "\n",
    "meas_dc_ext = meas.copy()\n",
    "simu_dc_ext = simu.copy()\n",
    "for n in names:\n",
    "    meas_dc_ext[n] = meas_dc_ext[n].extrapolate_to_dc(kind='linear')\n",
    "    simu_dc_ext[n] = simu_dc_ext[n].extrapolate_to_dc(kind='linear')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Microstripline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MSL_Thru_100` is a $L1$ long, $W1$ wide, $T$ thick copper microstripline on a $H$ height substrate with bottom ground plane.\n",
    "\n",
    "| Parameter | Value  |\n",
    "| :---      | :---   |\n",
    "| $L1$      | 100mm  |\n",
    "| $W1$      | 3.0mm  |\n",
    "| $T$       | 50umm  |\n",
    "| $H$       | 1.5mm |\n",
    "| Connector | Cinch, 142-0701-851|\n",
    "| Substrate | FR-4 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('tdr_measurement_vs_simulation/figures/MSL_100.jpg', width='50%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measurement vs simulation comparison\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(2,2,1)\n",
    "plt.title('Time')\n",
    "meas_dc_ext['P1-MSL_Thru_100-P2'].plot_z_time_step(0, 0, label='meas')\n",
    "simu_dc_ext['P1-MSL_Thru_100-P2'].plot_z_time_step(0, 0, label='simu')\n",
    "plt.xlim((-2, 3))\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.title('Frequency')\n",
    "meas_dc_ext['P1-MSL_Thru_100-P2'].plot_s_db(0, 0, label='meas')\n",
    "simu_dc_ext['P1-MSL_Thru_100-P2'].plot_s_db(0, 0, label='simu')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "z0 = 50\n",
    "t, ymea = meas_dc_ext['P1-MSL_Thru_100-P2'].s11.step_response(window='hamming', pad=0)\n",
    "ymea[ymea ==  1.] =  1. + 1e-12  # solve numerical singularity\n",
    "ymea[ymea == -1.] = -1. + 1e-12  # solve numerical singularity\n",
    "ymea = z0 * (1+ymea) / (1-ymea)\n",
    "t, ysim = simu_dc_ext['P1-MSL_Thru_100-P2'].s11.step_response(window='hamming', pad=0)\n",
    "ysim[ysim ==  1.] =  1. + 1e-12  # solve numerical singularity\n",
    "ysim[ysim == -1.] = -1. + 1e-12  # solve numerical singularity\n",
    "ysim = z0 * (1+ysim) / (1-ysim)\n",
    "plt.xlabel('Time (ns)')\n",
    "plt.ylabel('Relative error (%)')\n",
    "plt.plot(t*1e9, 100*(ysim-ymea)/ymea)\n",
    "plt.xlim((-2, 3))\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "delta = simu_dc_ext['P1-MSL_Thru_100-P2'].s_db[:,0,0] - meas_dc_ext['P1-MSL_Thru_100-P2'].s_db[:,0,0]\n",
    "f = simu_dc_ext['P1-MSL_Thru_100-P2'].f * 1e-9\n",
    "plt.xlabel('Frequency (GHz)')\n",
    "plt.ylabel('Delta (dB)')\n",
    "plt.plot(f, delta)\n",
    "plt.ylim((-20,20))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Surprisingly, the time domain results shows a very good agreement, within ±2%, while the frequency domain results exhibit a reasonable agreement only on the lower half of the frequencies.\n",
    "\n",
    "This is because the time domain base shape is mostly impacted by the low frequency data.\n",
    "\n",
    "There is a small offset between the time domain data, sign that the DC point is a bit different. The inductive peaks caused by the connector-to-microstripline transition are clearly visible and their simulation model agree well with measurements.\n",
    "\n",
    "To go further, it would be possible to measure up to a bigger frequency span, thus increasing the resolution and build the dut on a more frequency stable substrate, but this would be more costly. A cross section of the transmission line could also improve the knowledge of the actual geometry, including manufacturing tolerances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stepped impedance microstripline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MSL_Stepped_100` is a stepped microstripline made of $T$ thick copper on a $H$ height substrate with bottom ground plane. The section 1, 2, 3 and 4 from the left to the right, are described with length $Lx$ and width $Wx$, where $x$ is the section number.\n",
    "\n",
    "| Parameter | Value  |\n",
    "| :---      | :---   |\n",
    "| $L1$      | 50mm   |\n",
    "| $W1$      | 3.0mm  |\n",
    "| $L2$      | 20mm   |\n",
    "| $W2$      | 8.0mm (capacitive) |\n",
    "| $L3$      | 20mm   |\n",
    "| $W3$      | 1.0mm (inductive) |\n",
    "| $L4$      | 50mm   |\n",
    "| $W4$      | 3.0mm  |\n",
    "| $T$       | 50umm  |\n",
    "| $H$       | 1.5mm |\n",
    "| Connector | Cinch, 142-0701-851|\n",
    "| Substrate | FR-4 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('tdr_measurement_vs_simulation/figures/MSL_Stepped_140.jpg', width='75%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measurement vs simulation comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(2,2,1)\n",
    "plt.title('Time')\n",
    "meas_dc_ext['P1-MSL_Stepped_140-P2'].plot_z_time_step(0, 0, label='measurement')\n",
    "simu_dc_ext['P1-MSL_Stepped_140-P2'].plot_z_time_step(0, 0, label='simulation')\n",
    "plt.xlim((-1, 3))\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.title('Frequency')\n",
    "meas_dc_ext['P1-MSL_Stepped_140-P2'].plot_s_db(0, 0, label='measurement')\n",
    "simu_dc_ext['P1-MSL_Stepped_140-P2'].plot_s_db(0, 0, label='simulation')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "z0 = 50\n",
    "t, ymea = meas_dc_ext['P1-MSL_Stepped_140-P2'].s11.step_response(window='hamming', pad=0)\n",
    "ymea[ymea ==  1.] =  1. + 1e-12  # solve numerical singularity\n",
    "ymea[ymea == -1.] = -1. + 1e-12  # solve numerical singularity\n",
    "ymea = z0 * (1+ymea) / (1-ymea)\n",
    "t, ysim = simu_dc_ext['P1-MSL_Stepped_140-P2'].s11.step_response(window='hamming', pad=0)\n",
    "ysim[ysim ==  1.] =  1. + 1e-12  # solve numerical singularity\n",
    "ysim[ysim == -1.] = -1. + 1e-12  # solve numerical singularity\n",
    "ysim = z0 * (1+ysim) / (1-ysim)\n",
    "plt.xlabel('Time (ns)')\n",
    "plt.ylabel('Relative error (%)')\n",
    "plt.plot(t*1e9, 100*(ysim-ymea)/ymea)\n",
    "plt.xlim((-2, 3))\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "delta = simu_dc_ext['P1-MSL_Stepped_140-P2'].s_db[:,0,0] - meas_dc_ext['P1-MSL_Stepped_140-P2'].s_db[:,0,0]\n",
    "f = simu_dc_ext['P1-MSL_Stepped_140-P2'].f * 1e-9\n",
    "plt.xlabel('Frequency (GHz)')\n",
    "plt.ylabel('Delta (dB)')\n",
    "plt.plot(f, delta)\n",
    "plt.ylim((-10,10))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both time domain and the frequency domain results shows a reasonable agreement, within ±5% for time domain. The frequency domain results exhibit a good agreement, within ±dB only on the lower half of the frequencies, and the upper frequencies shows better agreement than in the previous [microstripline](#Microstripline) case.\n",
    "\n",
    "An explanation to the better frequency domain agreement is that the impedance steps induce such a discontinuities that more energi is reflected back, leading to a shorter effective length and an increased accuracy (more signal).\n",
    "\n",
    "The capacitive and inductive section of the stepped line are clearly visible and their simulation model agree well with measurements. The first connector effect is hidden by the plot scale, while the second is masked by the reflexions.\n",
    "\n",
    "To go further, it would again be possible to measure up to a bigger frequency span, thus increasing the resolution and build the dut on a more frequency stable substrate, but this would be more costly. A cross section of the transmission line could also improve the knowledge of the actual geometry, including manufacturing tolerances.\n",
    "\n",
    "Eventually, the stepped discontinuities could be made smaller, to produce a smaller effect on the overall measurement (each discontinuity has an influence on the following time domain signal shape)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dielectric effective relative permittivity and loss tangent characterisation based on multiline method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only two lines with different lengths are required for dielectric permittivity and loss tangent characterisation. Since we have measured the reflects too, we will use them instead of fake reflects. We don't have the switch terms either, but as we only extract the dielectric permittivity and loss tangent rather that doing a real calibration, there is no problem with that.\n",
    "\n",
    "Multiline calibration algorithm get rid of connectors effects by using multiple lengths of lines. The docstring explain that *At every frequency point there should be at least one line pair that has a phase difference that is not 0 degree or a multiple of 180 degree otherwise calibration equations are singular and accuracy is very poor.* These conditions will not be met with the chosen lines combination, but we will still be able to get a decent estimation of the dielectric permittivity and loss tangent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('tdr_measurement_vs_simulation/figures/MSL_100.jpg', width='50%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('tdr_measurement_vs_simulation/figures/MSL_200.jpg', width='100%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('tdr_measurement_vs_simulation/figures/MSL_Short_50.jpg', width='25%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Perform NISTMultilineTRL algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discard the *No switch term provided* warning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform NISTMultilineTRL algorithm\n",
    "line100mm = meas['P1-MSL_Thru_100-P2']\n",
    "line200mm = meas['P1-MSL_Thru_200-P2']\n",
    "short50mm = skrf.network.two_port_reflect(meas['P1-MSL_Short_50'], meas['P2-MSL_Short_50'])\n",
    "measured  = [line100mm, short50mm, line200mm]\n",
    "Grefls    = [-1]\n",
    "lengths   = [100e-3, 200e-3] # in meter\n",
    "offset    = [50e-3] # in meter\n",
    "cal = skrf.NISTMultilineTRL(measured, Grefls, lengths, er_est=4.5, refl_offset=offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Relative dielectric permittivity and loss tangent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `NISTMultilineTRL` calibration give the frequency-dependent effective relative dielectric of the geometry, which is a mix between air and substrate relative permittivity.\n",
    "\n",
    "Unfortunately, a single value of relative dielectric permittivity of the substrate at a given frequency is required for the microstripline media simulation, instead of the frequency-dependent effective value of the geometry we get from the calibration.\n",
    "\n",
    "To overcome this difficult situation, the microstripline model which will be used is fitted by optimization on the calibration results in such a way the difference between the frequency-dependent effective relative dielectric permittivity of both dataset is minimised.\n",
    "\n",
    "Additionally, a weighted contribution of the dielectric loss tangent is inserted in the optimization to minimise the difference between the measured modelled attenuation.\n",
    "\n",
    "The optimization results are dielectric $\\epsilon_r$ and $\\tan{\\delta}$ at 1 GHz and will be used in the dielectric dispersion model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# frequency axis\n",
    "freq  = line100mm.frequency\n",
    "f     = line100mm.frequency.f\n",
    "f_ghz = line100mm.frequency.f/1e9\n",
    "\n",
    "# the physical dimensions of the lines are known by design (neglecting manufacturing tolerances)\n",
    "W   = 3.00e-3\n",
    "H   = 1.55e-3\n",
    "T   = 50e-6\n",
    "L     = 0.1\n",
    "\n",
    "# calibration results to compare against\n",
    "ep_r_mea = cal.er_eff.real\n",
    "A_mea    = 20/log(10)*cal.gamma.real\n",
    "\n",
    "# starting values for the optimizer\n",
    "A     = 0.0\n",
    "f_A   = 1e9\n",
    "ep_r0 = 4.5\n",
    "tanD0 = 0.02\n",
    "f_epr_tand  = 1e9\n",
    "x0 = [ep_r0, tanD0]\n",
    "\n",
    "# function to be minimised\n",
    "def model(x, freq, ep_r_mea, A_mea, f_ep):\n",
    "    ep_r, tanD = x[0], x[1]\n",
    "    m = MLine(frequency=freq, Z0=50, w=W, h=H, t=T,\n",
    "            ep_r=ep_r, mu_r=1, rho=1.712e-8, tand=tanD, rough=0.15e-6,\n",
    "             f_low=1e3, f_high=1e12, f_epr_tand=f_ep,\n",
    "             diel='djordjevicsvensson', disp='kirschningjansen')\n",
    "    ep_r_mod = m.ep_reff_f.real\n",
    "    A_mod = m.alpha * 20/log(10)\n",
    "    return sum((ep_r_mod - ep_r_mea)**2)  + 0.1*sum((A_mod - A_mea)**2)\n",
    "\n",
    "# run optimizer\n",
    "res = minimize(model, x0, args=(freq, ep_r_mea, A_mea, f_epr_tand),\n",
    "               bounds=[(4.0, 5.0), (0.001, 0.1)])\n",
    "\n",
    "# get the results and print the results\n",
    "ep_r, tanD = res.x[0], res.x[1]\n",
    "print(f'epr={ep_r:.3f}, tand={tanD:.4f} at {f_epr_tand * 1e-9:.1f} GHz.')\n",
    "\n",
    "# build the corresponding media\n",
    "m = MLine(frequency=freq, Z0=50, w=W, h=H, t=T,\n",
    "         ep_r=ep_r, mu_r=1, rho=1.712e-8, tand=tanD, rough=0.15e-6,\n",
    "         f_low=1e3, f_high=1e12, f_epr_tand=f_epr_tand,\n",
    "         diel='djordjevicsvensson', disp='kirschningjansen')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calibration based values are plotted against modelled value as a sanity check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.subplot(2,2,1)\n",
    "plt.xlabel('Frequency [GHz]')\n",
    "plt.ylabel(r'$\\epsilon_{r,eff}$')\n",
    "plt.plot(f_ghz, ep_r_mea, label='measured')\n",
    "plt.plot(f_ghz, m.ep_reff_f.real, label='model')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.xlabel('Frequency [GHz]')\n",
    "plt.ylabel('A (dB/m)')\n",
    "plt.plot(f_ghz, A_mea, label='measured')\n",
    "A_mod = 20/log(10)*m.alpha\n",
    "plt.plot(f_ghz, A_mod, label='model')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.xlabel('Frequency [GHz]')\n",
    "plt.ylabel(r'$\\epsilon_{r,eff}$ error [%]')\n",
    "rel_err = 100 * ((ep_r_mea - m.ep_reff_f.real)/ep_r_mea)\n",
    "plt.plot(f_ghz, rel_err)\n",
    "plt.ylim((-2,2))\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.xlabel('Frequency [GHz]')\n",
    "plt.ylabel('$A$ error [%]')\n",
    "rel_err = 100 * ((A_mea - A_mod)/A_mea)\n",
    "plt.plot(f_ghz, rel_err)\n",
    "plt.ylim((-20,10))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The agreement between measurements and the model seems very reasonable. Relative error of $\\epsilon_{r,eff}$ stay within ±1% outside very low frequencies and relative error of $A$ is kept between ±10% on most of the range. Considering the shape of $A$, it is not possible to do much better with this model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connector effect characterization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `NISTMultilineTRL` calibration coefficients contain information about the connector characteristics, which is corrected by the calibration.\n",
    "\n",
    "These coefficients can be used to fit a connector model based on a transmission line section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Delay and attenuation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract connector characteristic from port 1 error coefficients\n",
    "conn = skrf.error_dict_2_network(cal.coefs, cal.frequency, is_reciprocal=True)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate connector delay with linear regression on unwrapped phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# connector delay estimation by linear regression on the unwrapped phase\n",
    "xlim = 9000 # used to avoid phase jump if any\n",
    "phi_conn = (np.angle(conn.s[:xlim,1,0]))\n",
    "z = np.polyfit(f[:xlim], phi_conn, 1)\n",
    "p = np.poly1d(z)\n",
    "delay_conn = -z[0]/(2*np.pi)\n",
    "print(f'Connector delay: {delay_conn * 1e12:.1f} ps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build connector model and compare it against calibration extracted data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc = DefinedAEpTandZ0(m.frequency, ep_r=1, tanD=0.02, Z0=50,\n",
    "                              f_low=1e3, f_high=1e18, f_ep=f_epr_tand, model='frequencyinvariant')\n",
    "\n",
    "Z0_conn = 50.0 # the actual connector characteristic impedance will be tuned later\n",
    "left = mc.line(delay_conn, 's', z0=Z0_conn)\n",
    "check = mc.thru() ** left ** mc.thru()\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "conn.plot_s_deg(1, 0, label='measured')\n",
    "check.plot_s_deg(1, 0, label='model')\n",
    "plt.ylabel('phase (rad)')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "conn.plot_s_db(1, 0, label='Measured')\n",
    "check.plot_s_db(1, 0, label='Model')\n",
    "plt.xlabel('Frequency (GHz)')\n",
    "plt.ylabel('Insertion Loss (dB)')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of connector model characteristics against calibration results shows a reasonable agreement. Calibration results exhibit some glitches that do not correspond to the expected physical behavior. They are caused by the calibration being close to singular due to the thru and line phase being a multiple of 180 degrees. Accuracy could be enhanced by feeding more distinct lines to the algorithm, but these are not manufactured yet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Characteristic impedance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have estimated connector delay and attenuation, but what about the characteristic impedance? This value is required to properly parametrize the transmission line section model.\n",
    "\n",
    "Optimization is used to find the characteristic impedance that minimize the difference between modelled and measured return loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s11_ref = conn.s[:,0,0]\n",
    "x0 = [Z0_conn]\n",
    "\n",
    "# function to be minimised\n",
    "def model2(x, mc, delay_conn, s11_ref):\n",
    "    Z0_mod = x[0]\n",
    "    conn_mod = mc.line(delay_conn, 's', z0=Z0_mod)\n",
    "    check = mc.thru(z0 = 50.) ** conn_mod ** mc.thru(z0 = 50.)\n",
    "    s11_mod = check.s[:,0,0]\n",
    "\n",
    "    return sum(absolute(s11_ref-s11_mod))\n",
    "\n",
    "# run optimizer\n",
    "res = minimize(model2, x0, args=(mc, delay_conn, s11_ref),\n",
    "               bounds=[(20, 100)])\n",
    "\n",
    "# get the results and print the results\n",
    "Z0_conn = res.x[0]\n",
    "print(f'Z0_conn={Z0_conn:.1f} ohm.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The modelled results are plotted against the calibration data, as a sanity check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn_mod = mc.line(delay_conn, 's', z0=Z0_conn)\n",
    "check = mc.thru(z0 = 50.) ** conn_mod ** mc.thru(z0 = 50.)\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "conn.plot_s_db(0,0, label = 'measured')\n",
    "check.plot_s_db(0,0, label = 'model')\n",
    "\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(check.f*1e-9, (check.s_db[:,0,0]-conn.s_db[:,0,0]))\n",
    "plt.ylabel('Delta (dB)')\n",
    "plt.xlabel('Frequency (GHz)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The delta in dB is quite big in low and high frequencies, but we will see when comparing measurement and simulation that the time domain reflectometry results are very decent with this value. The connector effects inductive peaks are well rendered in the case of the [microstripline](#Microstripline). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequency axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq = skrf.F(1,10e3,10000, 'mhz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Media sections with different geometries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 50 ohm segment\n",
    "MSL1 = MLine(frequency=freq, z0_port=50, w=W, h=H, t=T,\n",
    "        ep_r=ep_r, mu_r=1, rho=1.712e-8, tand=tanD, rough=0.15e-6,\n",
    "        f_low=1e3, f_high=1e12, f_epr_tand=f_epr_tand,\n",
    "        diel='djordjevicsvensson', disp='kirschningjansen')\n",
    "\n",
    "# Capacitive segment\n",
    "MSL2 = MLine(frequency=freq, z0_port=50, w=8.0e-3, h=H, t=T,\n",
    "        ep_r=ep_r, mu_r=1, rho=1.712e-8, tand=tanD, rough=0.15e-6,\n",
    "        f_low=1e3, f_high=1e12, f_epr_tand=f_epr_tand,\n",
    "        diel='djordjevicsvensson', disp='kirschningjansen')\n",
    "\n",
    "# Inductive segment\n",
    "MSL3 = MLine(frequency=freq, z0_port=50, w=1.0e-3, h=H, t=T,\n",
    "        ep_r=ep_r, mu_r=1, rho=1.712e-8, tand=tanD, rough=0.15e-6,\n",
    "        f_low=1e3, f_high=1e12, f_epr_tand=f_epr_tand,\n",
    "        diel='djordjevicsvensson', disp='kirschningjansen')\n",
    "\n",
    "# Connector transmission line media with guessed loss\n",
    "MCON = DefinedAEpTandZ0(m.frequency, z0_port=50, ep_r=1, tanD=0.025, z0=Z0_conn,\n",
    "        f_low=1e3, f_high=1e18, f_ep=f_epr_tand, model='frequencyinvariant')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulated devices under test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SMA connector\n",
    "conn_sma  = MCON.line(delay_conn, 's')\n",
    "\n",
    "\n",
    "# microstripline\n",
    "thru_simu = conn_sma ** MSL1.line(100e-3, 'm') ** conn_sma\n",
    "thru_simu.name = 'P1-MSL_Thru_100-P2'\n",
    "\n",
    "# stepped impedance microstripline\n",
    "step_simu = conn_sma \\\n",
    "           ** MSL1.line(50e-3, 'm') \\\n",
    "           ** MSL2.line(20e-3, 'm') \\\n",
    "           ** MSL3.line(20e-3, 'm') \\\n",
    "           ** MSL1.line(50e-3, 'm') \\\n",
    "           ** conn_sma\n",
    "step_simu.name = 'P1-MSL_Stepped_140-P2'\n",
    "\n",
    "# write simulated data to .snp files\n",
    "write_data = False\n",
    "if write_data:\n",
    "    step_simu.write_touchstone(dir='tdr_measurement_vs_simulation/simulation/')\n",
    "    thru_simu.write_touchstone(dir='tdr_measurement_vs_simulation/simulation/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After DC point extrapolation, and neglecting the effect of windowing, the time resolution of measurement is\n",
    "\\begin{equation*}\n",
    "Resolution = \\frac{1}{f_{span}}\n",
    "\\end{equation*}\n",
    "\n",
    "where $Resolution$ is the resolution in s and $f_{span}$ is the frequency span in Hz.\n",
    "\n",
    "In our case, with $f_{span} = 10 \\mskip3mu\\mathrm{[GHz]}$ \n",
    "\\begin{equation*}\n",
    "Resolution \\approx \\frac{1}{10^{10}} \\approx 100 \\quad \\mathrm{[ps]}\n",
    "\\end{equation*}\n",
    "\n",
    "With an effective relative dielectric permittivity of about 3.5, this give following physical resolution:\n",
    "\\begin{equation*}\n",
    "Resolution_{meter} = \\frac{Resolution \\dot{} c_0}{\\sqrt{\\epsilon{}_r}} \\approx \\frac{10^{-10} \\dot{} 3\\dot{}10^8}{\\sqrt{3.5} } \\approx 16 \\quad \\mathrm{[mm]}\n",
    "\\end{equation*}\n",
    "\n",
    "We will use reflexion measurement, so the actual range is divided by two, because the signal goes back and forth. Thus, the approximate distance resolution on our device under test will be 8mm. The discontinuities on the stepped impedance microstripline and the connector effects shall be well visible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measurement range limitation to avoid alias response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The measurement range should be set in such a way the true response of the device is seen without repetitions (alias).\n",
    "\n",
    "\\begin{equation*}\n",
    "Range = \\frac{c_0}{\\sqrt{\\epsilon{}_r} \\dot{} \\Delta{}_f}\n",
    "\\end{equation*}\n",
    "\n",
    "Where $Range$ is the range in m, $c_0$ the speed of light in m/s, $\\Delta{}_f$ the frequency step in Hz and $\\epsilon{}_r$ is the effective relative permittivity constant of the device under test.\n",
    "\n",
    "With the given measurement setup and considering $\\Delta{}_f = 1 \\mskip3mu\\mathrm{[MHz]}$, worst case $\\epsilon{}_r \\approx 5.0$, the range is\n",
    "\\begin{equation*}\n",
    "Range \\approx \\frac{3\\dot{}10^8}{\\sqrt{5} \\dot{} 10^6} \\approx 134 \\quad \\mathrm{[m]}\n",
    "\\end{equation*}\n",
    "Relative dielectric permittivity for FR-4 is approximately 5.0, the effective $\\epsilon{}_r$ of microstripline geometry will be smaller because of the air on top, which give a safety margin, the range being underevaluated.\n",
    "\n",
    "We will use reflexion measurement, so the actual range is divided by two, because the signal goes back and forth. However, our longest device being 200mm long, we can just quietly forget aliasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
